An efficient numerical method for calculating the entanglement of formation of 
arbitrary mixed quantum states of any dimension 
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It has been brought to our attention that this eprint duplicates earlier research by K. 
Audenaert et al: see quant-ph/0006128 or PRA 64 052304. We present a conjugate gradient 
method for calculating the entanglement of formation of arbitrary mixed quantum states of any 
dimension and with any bipartite division of the Hilbert space. The development of the gradient used 
by the algorithm, its implications for the number of states required in the optimal decomposition, 
and the way that conjugate gradient minimization has been adapted for this particular problem are 
outlined. We have found that the algorithm exhibits linear convergence for general mixed states, and 
that it correctly reproduces the known results for pairs of qubits and for isotropic states. The results 
of an example application of the code are discussed: calculating the entanglement of formation of a 
\^ + ) Bell state of two qutrits when one of those qutrits is subject to various decoherence channels. 
The results for qutrits are contrasted with those for qubits: for the types of decoherence considered 
here, qutrit entanglement appears to be more robust than qubit entanglement. 

PACS numbers: 03.67.-a,03.65.Ud,05.30.-d 



Introduction. The entanglement of formation Ep 
characterizes the amount of entanglement present in a 
mixed state p. The definition of Ep is known: it is the 
average entanglement of the pure states {|Vi)} m the 
decomposition of p minimized over all possible decom- 
positions {pi, IV'i)} °f P- However, except for the special 
case of two qubits [|, no formula for Ep has been found, 
nor a prescription for obtaining the decomposition which 
realizes the minimum. 

In this Letter, we describe the development of a gra- 
dient for the average entanglement E av , and a conjugate 
gradient algorithm that uses it to calculate Ep numeri- 
cally. The algorithm successfully minimizes the average 
entanglement of random two qutrit (3 x 3) mixed states 
within a few hours on an ordinary desktop computer. 
The algorithm has been tested against the known result 
for general states of two qubits, and against exact results 
for Ep for the special class of 'isotropic' mixed states of 
two qutrits. Although the examples discussed in this pa- 
per involve a division of the Hilbert space into subspaces 
of equal dimension, the algorithm works equally well for 
any division of the Hilbert space. For example, it could 
be used to study the entanglement between a qutrit and 
quqit (d=4 system) whose joint state is mixed. 

The entanglement of formation. A particular pure 
state decomposition of a mixed state p shared between 
subsystems A and B can be written 



P = 



= ^Pik = ^Pi\$i)($%\=^\i>i)$i\ (1) 
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where the {IV'i)} are subnormalized states defined by 

IV'i) := y/pl\ipi) 0- The average entanglement of the 
decomposition, in ebits, is 

i 

= -log 2 eJ2Pi^A[ptlnpf} (2) 



where Ei is the entanglement of IV'i)- For convenience, 
the derivation works with natural logs. The entangle- 
ment of formation is defined as the minimum of E av over 
all possible decompositions of p: 



Ep(p) := min E av 



(3) 



The decomposition which gives the minimum is described 
as 'optimal'. It is possible to make unitary transforma- 
tions between different decompositions of a given mixed 
state 0: 



iVi) 



(4) 



Although these decompositions all generate the same 
mixed state, they have different values of the average 
entanglement. 

A gradient for the average entanglement. Now write 



U = exp(ie9) 



(5) 



where e is a real parameter, and 9 is a Hermitian matrix. 
All derivatives are evaluated at e — 0. The gradient of 
pf is easily obtained: 



dp{ 



1 d 



Pi de 



= -x (Tr B (| ViKViD) - 3 ^Tr B (|Vi)(Vil06) 
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The gradient of Tr^l/i^ In pf] requires a little more 
work. If a labels an eigenvalue of pf , and \oti) is the 
corresponding eigenket, 



= Mln^f 



(7) 



where we have made use of the invariance of Tr(pf) and 
the Hellman-Feynman theorem 0- C3- 01 ■ Now we can 
calculate the promised expression for the gradient: 



dE a , 



-log 2 e W^Tr^lnpf] + p 4 Tr[ln p A ^ 



dc 



- lo g2 e X! Tr ' 4 



lnpf-Tr B |^){^ 



(8) 



The result can be further simplified by calculating the 
gradient of Pipf — Tr B (|-0i) and making use of equa- 
tion 



U 



dc 



(Pipi) 



Je9 



dc 



dc 



U, 



e=0 



d 

~de 



dE„ 



dc 



; -iTr B |^)(^|) 



-iE% Tr ^ 
E 



log 2 p? - \og 2 pf Tr B |^-)(^| 



(10) 



where g is the matrix of gradient elements with respect 
to the space of generators {9}. 

Number of states required in the optimal decomposi- 
tion. Suppose we add to our set of states a new one, 
sub-normalized to zero, and allow unitary transforma- 
tions among this new, enlarged, set. Such transforma- 
tions automatically leave the density matrix invariant, 
but it seems possible that they might reduce the average 
entanglement. However, equation (|10[l shows that they 
do not, at least to first order, since the trace over B on the 
RHS will be zero if either or \ipj) is sub-normalized 
to zero. Therefore no 'padding' with zero-norm states 



is required to produce the correct result. This strongly 
suggests that the number of states needed is equal to the 
rank of the density matrix p, instead of the previously 
accepted bound of the square of the dimension of p. 

Existence of local and global minima In general, a 
conjugate gradient procedure will converge to a local, 
rather than a global, minimum. However, Prager 
has shown that any local minimum of Ep is also a global 
minimum. We have encountered points (presumably lo- 
cal maxima or inflexion points) other than the optimal 
decomposition at which the gradient is zero: when this 
occurs we find we can restart the conjugate gradient al- 
gorithm after a small number of random moves. 

Using the gradient to minimize the average entangle- 
ment. The code described here uses the conjugate gra- 
dient algorithm to minimize E av with respect to the space 
of all possible unitary transformations between some ini- 
tial decomposition of p (the choice of which is entirely 
arbitrary) and the current decomposition. The initial de- 
composition is chosen without prejudice to be the 'eigen- 
state' decomposition of p, ie. consisting of pure states 
which are the eigenvectors of p and probabilities which 
are the eigenvalues. We write U = exp(iH) and move 
through the space of Hermitian matrices {H} by con- 
structing conjugate gradient moves from the gradient in- 
formation provided by g. The initial unitary transforma- 
tion is the identity I, whose corresponding H is the null 
matrix 0. 

The standard formulation of the conjugate gradient al- 
gorithm [t| operates upon vectors. Therefore the point 
in H-space used is the 'flattened' vector consisting of the 
unique elements of H. The gradient matrix is flattened in 
the same way. Conjugate directions are chosen according 
to the Fletcher- Reeves- Polak-Ribiere algorithm [8] . The 
gradient information is also made use of by the modi- 
fied version of Brent's method [9j that performs the line 
minimizations. The end result of the algorithm is the 
unitary transformation which takes us from the 'eigen- 
state' decomposition to the optimal decomposition. An 
heuristic test for convergence is performed after the al- 
gorithm terminates by running Monte Carlo against the 
final decomposition, using an exponentially wide range 
of step sizes in random directions through unitary trans- 
formation space: we fail to find further steps down to 
within the target precision. This does not prove that the 
optimal decomposition has been obtained, but is strongly 
suggestive of that conclusion. Note that the gradient for 
E av is constructed with respect to small unitary transfor- 
mations of the current decomposition, and therefore the 
unitarity of U is preserved to machine precision through- 
out the execution of the conjugate gradient algorithm. 

Performance against two-qubit mixed states The al- 
gorithm's estimate of Ep converges with the Wootters 
Ep to 14 decimal places typically in less than 100 itera- 
tions. 

Performance against two-qudit isotropic mixed states 
An isotropic mixed state is a convex mixture of the maxi- 
mally entangled Bell state and the maximally mixed 
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state pi = I: 

p ■= ^-r(i - i* + x* + d + fi*+><*+i (ii) 

Unlike for more general mixed states, a formula for Ep 
exists for isotropic states of two qudits |l(J . These states 
are thus well-suited as a check for the results of the algo- 
rithm. In the case of two-qubit isotropic states, a partic- 
ular difficulty arises which is that E &v generally has a sta- 
tionary point at the eigenvector decomposition. It is thus 
necessary to move the initial decomposition away from 
this stationary point with a few random Monte Carlo 
moves, after which the conjugate gradient algorithm suc- 
cessfully obtains the correct minimum. For two-qutrit 
isotropic states we observe perfect correspondence be- 
tween the results of the algorithm and the value given by 
the formula. 

Performance against two-qutrit random mixed states. 
The convergence behaviour when minimizing a sample 
random general mixed state of two qutrits (ie. dim(7i) = 
9) is shown in Figure It can be seen that convergence 
is linear, with approximately 600 iterations required per 
significant figure. We have found a few examples of states 
with highly degenerate eigenvalue spectra (e.g. isotropic 
states) for which convergence may be slower than linear. 

Modelling locally depolarized Bell states of two qutrits. 
An example of a previously intractable yet conceptually 
simple problem is local depolarization of one qudit in- 
volved in an entangled state. Consider a \^ + ) state of 
two systems of dimension d (qudits), defined by 

l* + > : =-73l>>- ( 12 ) 

If one qudit - say qudit A - decoheres with a certain 
probability p, the entanglement of the overall state is 
reduced in a non-trivial manner. This corresponds to a 
physical situation in which one party (Bob) prepares two 
qudits in an entangled state, then passes one qudit to 



Alice via a noisy channel. How much entanglement do 
Alice and Bob now share? 

In this exercise, three possible types of decoherence are 
considered: the bitflip channel, the depolarizing channel, 
and both channels in succession with the same probabil- 
ity. (As the operations commute in this case, the order 
of bitflip and depolarization operations in the last case is 
irrelevant.) The results for a |\E' + ) state of two qubits can 
be readily calculated using the Wootters formula. They 
are shown in Figure [21 For the depolarizing channel, the 
results are consistent with the well-known separability 
of mixed states in the neighbourhood of the maximally 
mixed state [llj . 

To perform the same exercise for a |* + ) state of two 
qutrits (d = 3), the bitflip channel is extended to qutrits 
by using an operator sum representation consisting of 
\/l _ pl> vW2X up , ^Jp~j2~K d own-, where p represents the 
probability of a flip and the latter two operators respec- 
tively raise and lower a basis state. The depolarizing 
channel is extended to qutrits using the operator-sum 
representation in |l2j with d = 3. The results of apply- 
ing these operations (tensored with the identity on qutrit 
B, which does not decohere) are shown in Figure |3J The 
results are plotted in etrits (I etrit = log 2 (3) ebits) to 
facilitate comparison with the results for qubits. (Note 
that for the depolarization-only case, the results can also 
be obtained using the formula for Ep for isotropic states 
- we have found they are the same). They are similar 
to those for qubits, however it can be seen that for all 
types of channel, the proportionate reduction in entan- 
glement at a given probability is lower for qutrits. It thus 
appears that qutrit entanglement is more robust with re- 
spect to these types of decoherence than qubit entangle- 
ment. This conclusion has potential significance for the 
design of quantum computers and communication chan- 
nels. 
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FIG. 1: Convergence of min(i? av ) for a random two-qutrit mixed state using conjugate gradient algorithm. 
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FIG. 2: Entanglement of formation (ebits) of \^ + ) of two qubits vs probability of qubit A depolarizing through various channels 
(calculated using Wootters formula). 
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FIG. 3: Entanglement of formation (etrits) of |^ + ) of two qutrits vs probability of qutrit A depolarizing through various 
channels (calculated using conjugate gradient algorithm) . 



Qutrit A [lipped - 
Qutrit A depolarized 
Qutrit A Hipped then depolarized - 




p (decoherence event/s on qutrit A) 



